# TODO: Add comment
# 
# Author: guochun
###############################################################################

rpp=function(N,plotdim,rthin=FALSE,f=NULL){
	n=round(sqrt(N/2),0)
	x=seq(0,plotdim[1],length.out=n)
	y=seq(0,plotdim[2],length.out=n)
	
	nn=length(y)
	x=rep(x,nn)
	y=rep(y,each=nn)
	if(rthin){
		p1=ppp(x=x,y=y,window=owin(c(0,plotdim[1]),c(0,plotdim[2])))
		p1=rthin(p1,f)
		x=p1$x
		y=p1$y
	}
	x=x+rnorm(length(x),0,1)
	y=y+rnorm(length(y),0,1)
	xx=x+rnorm(length(x),0,1)
	yy=y+rnorm(length(y),0,1)
	mark1=runif(length(x),3,5)
	mark2=mark1-2
	re=ppp(x=c(x,xx),y=c(y,yy),window=owin(c(0,plotdim[1]),c(0,plotdim[2])),marks=c(mark1,mark2))
	return(re)
}

rmpp=function(pp,r,corfun,markrange=c(1,10)){
	corfun=eval(corfun)
	#pp=rpoispp(lambda,win=window)
	N=pp$n
	marks=rep(NA,times=N)
	for(i in 1:N){
		centx=pp$x[i]
		centy=pp$y[i]
		neibs=pp$x>=(centx-r) & pp$x<=(centx+r) &
	          pp$y>=(centy-r) & pp$y<=(centy+r)
	    if(any(neibs)){
			selmarks=marks[which(neibs)]
			if(any(!is.na(selmarks))){
				selmarks=selmarks[!is.na(selmarks)]
				if(length(selmarks)==1){
					marks[i]=selmarks+markrange[2]/4
					if(marks[i]>markrange[2])
						marks[i]=marks[i]-markrange[2]
				}else{
					p=density(selmarks,from=markrange[1],to=markrange[2])$y
					p=max(p)-p
					#mm=mean(selmarks,na.rm=T)
					#p=((mds-mm)^2)^(1/3)
					#p=p/sum(p)
					marks[i]=sample(seq(markrange[1],markrange[2],length.out=length(p)),1,prob=p)
				}
				
			}else
				marks[i]=runif(1,markrange[1],markrange[2])
		}else{
			marks[i]=runif(1,markrange[1],markrange[2])
		}
	}
	marks(pp)=marks
	return(pp)
}

envelopeKm=function(ppp,ppp2=NULL,mf,r=NULL,plot=FALSE,...){
	re=Km(ppp,mf,r=r)
	re.g=pcfConvert(re)
	
	temp=numeric()
	rep=39
	X=ppp
	if(is.null(ppp2)){
		for(i in 1:rep){
			#X$x=runif(X$n,0,1000)
			#X$y=runif(X$n,0,500)
			X$marks=sample(X$marks)
			re.k=pcfConvert(Km(X,mf,r=r))
			temp=c(temp,re.k[,2])
		}
	}else{
		re2=Km(ppp2,mf,r=r)
		re.g2=pcfConvert(re2)
		
		Nchange=ppp2$n-ppp$n
		for(i in 1:rep){
			del=sample(1:ppp2$n,Nchange)
			X2=ppp2[-del]
			re.k=pcfConvert(Km(X2,mf,r=r))
			temp=c(temp,re.k[,2])
		}
	}
	
	dim(temp)=c(dim(re.k)[1],rep)
	re.max=apply(temp,1,max)
	re.min=apply(temp,1,min)
	re.g$max=re.max
	re.g$min=re.min
	
	if(plot){
		x=re.g[,1]
		plot(x=range(x),y=range(c(re.g[,2],re.g[,3],re.g[,4])),type="n",...)
		lines(x=x,y=re.g[,2])
		lines(x=x,y=re.g[,3],col=2,lty=2)
		lines(x=x,y=re.g[,4],col=2,lty=2)
		if(!is.null(ppp2)){
			lines(x=x,y=re.g2[,2],col=2)
		}
	}
	if(is.null(ppp2))
	    return(re.g)
	else
		return(list(re.g,re.g2))
}